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The interest in studying metabolic alterations in cancer and their potential role as novel targets for 
therapy has been rejuvenated in recent years. Here, we report the development of the first genome- 
scale network model of cancer metabolism, validated by correctly identifying genes essential for 
cellular proliferation in cancer cell lines. The model predicts 52 cytostatic drug targets, of which 
40% are targeted by known, approved or experimental anticancer drugs, and the rest are new. It 
further predicts combinations of synthetic lethal drug targets, whose synergy is validated using 
available drug efficacy and gene expression measurements across the NCI-60 cancer cell line 
collection. Finally, potential selective treatments for specific cancers that depend on cancer 
type-specific downregulation of gene expression and somatic mutations are compiled. 
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Introduction 

The interest in studying cancer metabolism has recently grown 
in light of the decreasing number of newly released anticancer 
drugs, the development of metabolite profiling technologies 
and metabohc network databases, and due to increased efforts 
to target metabohc diseases by the pharmaceutical industry. 
During tumor development, cancer cells modify their meta- 
bohsm to meet the requirements of cellular prohferation, thus 
facilitating the uptake and conversion of nutrients into 
biomass (Vander Heiden et aZ, 2009} . Numerous studies have 
shown similar metabohc alterations occurring across tumor 
cells, including changes in glucose metabohsm that give rise to 
the Warburg effect, and an increase in biosynthetic activities 
(such as nucleotide, lipids and amino-acid synthesis), im- 
portant for cellular prohferation (DeBerardinis et aU 2008; 
Tennant et aU 2009; Vander Heiden et aU 2009} The latter 
are regulated by several signaling pathways implicated in 
ceh prohferation, in which cancer-associated mutations lead 
to accelerated growth (Christofk et aU 2008). The similarity 
in metabohc activity across cancer types is evident in 
experimental measurements of metabolic flux and from 
analyzing multiple high-throughput molecular data sources, 
including gene expression, metabolomics and phenotypic 
assays of ceh line growth fohowing shRNA-facilitated gene 
silencing (knockdown) (Supplementary Information). Meta- 
bohc commonalities across tumors have important practical 
imphcations, where drugs, such as antimetabolites that target 
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nucleotide biosynthesis are used to fight a variety of different 
malignancies. 

The observation that many types of cancer cehs adapt their 
metabohsm to facilitate biomass formation to enable proh- 
feration suggests that it may be possible to predict characte- 
ristic alterations in cancer metabohsm via genome-scale 
computational modeling approaches that have been success- 
fully used in the past to predict the metabohc state of fast- 
growing microorganisms (Price et aU 2004). In particular, flux 
balance analysis (FBA) is a constraint-based modeling (CBM) 
approach that is suitable for modeling cancer metabohsm as it 
assumes that a ceh is under selective pressure to increase its 
growth rate, and hence searches for metabohc flux distribu- 
tions that produce essential biomass precursors at high rates 
(Price et aU 2003). A fundamental step toward large-scale 
modeling of human metabohsm has been taken in recent 
studies (Duarte et aU 2007; Ma et aU 2007), which recon- 
structed a generic (non-tissue specific) human metabohc 
network based on an extensive evaluation of genomic and 
bibhomic data. The potential clinical utility of a human 
metabohc model has already been demonstrated, showing 
its ability to identify reactions causally related to hemolytic 
anemia (Duarte et aU 2007), to predict potential drug targets 
for treating hypercholesterolemia (Duarte et aU 2007), disease 
co-morbidity (Lee et aU 2008) and tissue specificity of disease 
genes (Shiomi et aU 2008), in identifying diagnostic bio- 
markers for inborn errors of metabohsm (Shiomi et aU 2008). 
A recent study by Li et al (2010) have further utilized the 
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human metabolic model as a basis for the prediction of novel 
targets for known anticancer drugs (Li et al, 2010} . 

Results 

A genome-scale model of cancer metabolism 
correctly predicts cancer growth-supporting 
genes 

Here, we construct the first large-scale model of cancer 
metabohsm that aims to capture the main metabohc functions 
common across many cancer types. To this end, we integrate 
the human metabohc model of Duarte et al (2007) with 
cancer gene expression data, utilizing a variant of our recent 
computational method for the automatic reconstruction of 
human tissue metabohc models, termed MBA (Model Building 
Algorithm) (Jerby et aZ, 2010] (see Materials and methods). 
In brief, we begin by assembling an initial core set of 197 
metabohc enzyme-coding genes that are highly expressed 
across 90% of all cancer cell lines in the NCI-60 collection 
(Grever et aU 1992). Then, applying MBA, a minimal set 
of additional reactions from the human metabolic model that 
are needed to activate the reactions associated with this initial 
core set is added, obtaining a cancer metabohsm model that 
is consistent, that is, each of the core reactions can poten- 
tially carry a non-zero metabolic flux, within a global flux 
distribution satisfying stoichiometric, mass balance and reac- 
tion directionality constraints. The minimal set of added reac- 
tions was chosen such that it enables the steady-state 
production of biomass compounds required for cellular proh- 
feration (based on prior knowledge of their relative concentra- 
tions; see Materials and methods). This model reconstruction 
approach, accounting for both a core set of highly expressed 
genes along with a cellular proliferation constraint extends 
upon previous model reconstruction methods (Shlomi et al, 
2008; Jerby et al, 2010). The resulting cancer metabohc model 
includes 772 reactions and 683 genes (see Materials and 
methods; Supplementary Table S8). Applying FBA to this 
model enables the prediction of the metabohc state of cancer 
cells (their prohferation rate and flux distribution) across 
different gene knockdowns and modeling the effects of drug 
apphcations on a large scale. 

To vahdate the cancer network model generated, we 
analyzed it via FBA to predict growth-supporting genes whose 
knockdown would significantly reduce cellular proliferation 
rate (see Materials and methods), that is, we knocked down all 
genes one at a time, and recorded the resulting simulated 
growth rate in the model. Overall, we obtain a set of 199 
growth-supporting genes, which are reassuringly ranked as 
highly essential based on shRNA gene silencing data (Luo et al, 
2008) (Kolmogorov-Smirnov (KS) P-value=0.0045; Supple- 
mentary Figure SI) . The reference shRNA data set consists of a 
hst of genes, ranked according to the survival rate of 12 cancer 
cell lines after these genes are knocked down, thus denoting 
the genes' experimentally measured contribution to cancer 
growth. As a further validation, we compared the predicted set 
of growth-supporting genes with a set of genes in which cancer 
somatic mutations have been reported in the literature (Futreal 
et al, 2004), finding a high level of enrichment (hyper- 
geometric P-value of 0.002). For comparison, we evaluated 
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the predictive performance of the gene expression data and the 
human metabohc model separately, finding that either source 
is a worse predictor of both shRNA gene essentiality (P-values 
are 0.075 and 0.029, respectively) and cancer mutations 
(P-values are 0.0077 and 0.025, respectively) compared with 
the generic cancer model (Materials and methods; Supple- 
mentary Information) . 

While we described above the reconstruction of a generic 
cancer metabohc model, capturing major metabohc altera- 
tions that are common across cancer types, a similar compu- 
tational approach can be utilized to develop more fine-tuned 
metabohc models of specific cancer types. To demonstrate this 
potential, we apphed our cancer model building approach to 
reconstruct a model of non-smah ceh lung cancer metabohsm 
utilizing multiple gene expression data sets (Supplementary 
Information). The growth-supporting genes predicted by 
the lung cancer model are ranked as highly essential based 
on shRNA gene silencing data measured in this ceh line 
(KS P-value=0.025), outperforming the predictions made by 
the generic cancer model for this ceh line (Supplementary 
Information) . 

Predicting cytostatic anticancer targets 

To identify which of the above-identified growth-supporting 
genes may be considered viable anticancer targets, we further 
aimed to predict whether their knockdown is expected to be 
toxic to non-dividing cehs or damage the proliferation of 
normal cehs. To simulate the effect of these gene knockdowns 
on normal, non-dividing cehs, we predicted the potential 
damage to ATP production (utilizing the entire human 
metabohc network), a vital biochemical function that must 
be preserved in every ceh. Based on these predictions, we 
define a cytostatic score for each gene, representing the extent 
to which its knockdown reduces cancer growth compared with 
its effect on ATP production in healthy cehs (see Materials and 
methods) (with a cytostatic score of 1 denoting a non-toxic 
target that completely eliminates cancer growth without 
affecting ATP production in healthy cehs). The resulting 
distribution of cytostatic scores is bimodal; out of the 199 genes 
that are predicted to be growth supporting in the cancer model, 
52 have a high cytostatic score (above 0.9), and the remaining 
147 genes have a low score (below 0.6; Figure lA; Supple- 
mentary Table S2; Supplementary Information). As an addi- 
tional method to predict the effect of these knockdowns on the 
metabohsm of healthy, non-dividing cehs, we tested the effects 
of knockdowns in a model of liver metabohsm (Jerby et al, 
2010) (specifically on normal urea secretion, glycogenesis, 
glycogenolysis, gluconeogenesis and bilirubin uptake), ruling 
out one of these drug targets (CMPKl) as potentially damaging 
normal liver uptake of bilirubin (see Materials and methods; 
Supplementary Table S2) . Notably, in the absence of detailed 
metabohc networks for a wide range of different human 
tissues, it is currently impossible to strictly rule out that the 
predicted targets would not damage metabolic functions of 
other healthy tissues. As a final screening step, we studied how 
gene knockdowns would affect prohferation of healthy cehs. 
To this end, we apphed a standard FBA analysis on the entire 
human metabolic network model, aiming to identify growth- 
supporting genes as done for the cancer model. We found that 
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Cytostatic score Pathway 

Figure 1 Cancer selectivity and pathway association of predicted growth-supporting genes. (A) Distribution of selectivity scores for the set of 199 predicted growth- 
supporting genes. (B) Pathway association of the highly cytostatic growth-supporting genes (cytostatic score > 0.9), showing for each pathway the number of predicted 
genes that are known targets of current anticancer drugs, the number of known targets of drugs that are currently used for non-cancer indications and entirely novel gene 
targets, that is, genes without any currently known drugs that target them. For each pathway, the number of missed predictions, that is, known anticancer drug targets 
that are not predicted to be highly selective, is also shown. 



the knockdown of 49 out of the 52 high cytostatic scoring 
targets is likely to also damage prohferation of normal cells 
(Supplementary Table S2}, suggesting that the targeting of 
these genes would cause similar side effects to current 
cytostatic drugs (Partridge et aU 2001). 

The 52 targets with high cancer cytostatic scores contain 
8 out of 24 known targets of the 14 FDA-approved metabohc 
anticancer drugs found in DrugBank (Wishart et al, 2008} 
(Supplementary Tables SI, S2 and S5}, representing a highly 
significant enrichment (hypergeometric P- value <6xl0~''}, 
testifying to predictive power of the model. Interestingly, 
the model misses the prediction of the knockdown effects 
on prohferation of 16 known drug targets. These include 
four genes involved in the metabohsm of signaling molecules 
(e.g. aminoglutethimide, which targets a hormone-signaling 
mechanism), whose effect on cellular growth is outside the 
scope of the model; seven additional genes that are targeted by 
antimetabolites (Mercaptopurine, Thioguanine and Capecita- 
bine) having multiple metabolic targets (Supplementary 
Table S5} . Finally, four additional genes are correctly predicted 
as part of synthetic lethal simulations described below, 
such that overall, out of all the metabolically active drug 
targets currently known, only one missed prediction remains 
unaccounted for. Moving from the gene to the drug level 
and simulating drug treatments by concurrently knocking 
down multiple genes targeted by the same drug, we success- 
fully identify all three antimetabolites as having a high 
cytostatic score, correctly predicting 11 (78%) out of all 
14 FDA- approved metabolic anticancer drugs (Supplementary 
Table S5}. 

Interestingly, the set of 52 predicted highly cytostatic genes 
contains 13 additional genes that are targeted by existing non- 
cancer drugs. Remarkably, all of these predicted drugs but one 
are currently under experimental testing for cancer therapy 
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(Supplementary Table S2}: eight of these genes are targeted by 
drugs already approved by the FDA for the treatment of various 
non-cancer diseases (e.g. antimalarial treatment, hyper- 
cholesterolemia, obesity, etc.) and five genes are targeted by 
experimental drugs. The remaining 31 (out of the original 52] 
cytostatic genes are currently without any known drug that 
currently targets them, and form interesting candidates for 
potential drug targeting. Notably, these novel drug targets 
span a much broader set of potential anticancer target path- 
ways; yet, reassuringly many of them participate in the same 
metabolic pathways as the known and experimental anti- 
cancer drug targets (Figure IB; Supplementary Table S2}. 
For example, we predict five highly selective drug targets in 
sphingolipids metabohsm (a subclass of fatty acid metabohsm; 
currently not targeted by any FDA-approved drugs), whose 
targeting was recently shown to represent a promising 
new approach to treat pancreatic cancers (Guillermet-Guibert 
et al 2009). 

Predicting synthetic lethal gene targets 

The concept of synthetic lethality promises to have a central 
role in the future selection of anticancer drug targets (Kaelin, 
2005). Synthetic lethal genes are common in metabohc 
networks due to their highly robust nature (resulting from 
backup mechanisms such as isozymes and alternative path- 
ways; Stelling et aU 2002), and their identification via FBA was 
successfully demonstrated in microbial networks (Deutscher 
et aU 2006; Harrison et aU 2007] . To study synthetic lethal drug 
targets, we systematically simulated all double gene knock- 
downs in the cancer model. We assigned each gene pair with 
a synergy score, reflecting the additional drop in prolife- 
ration rate below the minimal rate achieved by its indivi- 
dual single knockdowns (in accordance with the commonly 
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used 'highest single agent' synergy model; Berenbaum, 1989). 
The analysis resulted in 342 synthetic lethal gene pairs with a 
synergy score >0.5. 

As a vahdation for the plausibility of the predicted 
synergistic gene pairs, we find that they are significantly 
enriched with genetic interactions between the corresponding 
yeast orthologs (hypergeometricP-value=0.028; see Materials 
and methods), considering a large-scale data set of genetic 
interactions in yeast (Costanzo et aU 2010). Notably, the 
significance of the enrichment is quite remarkable, consi- 
dering that the conservation level of genetic interactions 
between distant species is considered to be rather low (as 
demonstrated, e.g., between Saccharomyces cerevisiae and 
Caenorhabditis elegans; Tischler et aU 2008). As a second 
validation, we investigated whether drugs that target a single 
gene (participating in a predicted synthetic lethal pair) 
indeed have higher efficacy in cell lines in which the 
synergistic gene is lowly expressed. To this end, we analyzed 
together growth inhibition measurements for 11 metabohc 
drugs for which such data exists (Lee et aU 2007} and 
gene expression measurements over all NCI-60 cell lines 
(Grever et al 1992} (Figure 2A}. We find that drugs that 
target one of the genes in a predicted pair indeed exert a 
significantly stronger inhibitory effect on prohferation in 
the cell lines in which their paired target gene has a low 
expression level (P-value=0.02}, supporting our predictions 
(Figure 2B}. 
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Figure 2 The relation between drugs' efficacy and the expression of their 
paired target genes. (A) A schematic illustration of the expected relation between 
an efficacy of a drug and the expression of genes that have a synthetic lethal 
interaction with the drug's target. (B) The distribution of efficacy-expression 
correlation values for the set of predicted synthetic lethal drug-gene pairs, and for 
a background distribution of random pairs. The predicted synthetic lethal pairs 
show a significant anticorrelation compared with the random pairs (Wilcoxon 
P-value=0.02). 



The targeting of both synthetic lethal genes via 
combination therapy 

To identify which of the synthetic lethal genes may be further 
considered for combinatorial drug therapies, we further 
predicted whether their joint knockdown is expected to be 
toxic for non-dividing cells or damage the prohferation 
of normal cells (as performed above for the single targets}. 
We find that 133 gene pairs have high cytostatic score 
(>0.6; computed when they are both knocked down; 
Figure 3 A; see Materials and methods}, suggesting (within 
the screening limitations pointed above} that their double 
knockdown may not be toxic to non-dividing cells. Exploration 
of the metabohc pathways in which the predicted highly 
cytostatic gene pairs participate (Figure 3B; Supplementary 
Table S4} reveals that both the pentose-phosphate path- 
way (PPP} and pyruvate metabolism have a central role, with 
45% of the gene pairs having at least one target in these 
pathways. PPP is involved in the production of both NADPH 
(necessary for lipid and nucleotide biosynthesis as well 
as protection against oxidative stress} and ribose 5-phosphate 
(a nucleotide precursor}, making it an attractive target for 
cancer therapy, though currently there are still no inhibitors 
for this pathway in clinical trials. In pyruvate metabo- 
hsm, pyruvate carboxylase (PC; participating in 13 pairs} 
is notable; it is involved in one of the two major pathways 
that contribute to anaplerosis (i.e. the replenishment of 
TCA cycle intermediates}, the other pathway being glutami- 
nolysis. Indeed, several previous attempts have already 
been made to target glutaminolysis (Rosenfeld and Roberts, 
1981}, further supporting the plausibility of targeting PC 
as well. 

Predicting how the knockdown of the synthetic lethal targets 
would affect prohferation of normal cells, we find that in a 
marked contrast to the single target analysis, the knockdown 
of 99 of the pairs is not expected to adversely damage growth 
of normal cells (Supplementary Table S3}. The fact that 
synthetic lethal genes are predicted to achieve selective 
targeting of cancer cells is in accordance with previous 
experimental findings (Lehar et al, 2009}. For example, 
we predict that the double knockdown of glycine hydroxy- 
methyltransf erase (SHMTl} and alanine-glyoxylate trans- 
aminase (AGXT} to differentially inhibit proliferation of cancer 
versus normal cells. While these knockdowns would block 
glycine biosynthesis in cancer cells (from serine and glyoxy- 
late, respectively; Supplementary Figure S2A}, their knock- 
down in normal cells is predicted to spare glycine production 
through an alternative pathway involving choline degrada- 
tion. The latter pathway was not included in the cancer 
model due to low expression level of its member genes (while 
choline utilization in the cancer model proceeds via choline 
kinase toward phosphocholine biosynthesis, in agreement 
with experimental data; Nakagami et aU 1999}. Other 
examples include synthetic lethal genes that are predicted to 
differentially inhibit nutrient uptake rates in cancer versus 
normal cells, due to the exclusion of several lowly expressed 
membrane transporters in the cancer model (Supplementary 
Figure S2B and C}. 

Predicting the effect of the 133 double gene knockdowns on 
liver metabohsm resulted in eight pairs whose targeting could 
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Figure 3 Cytostatic scores, synergy and pathway association of predicted target gene pairs. (A) The distribution of predicted drug target pairs across their synergy and 
selectivity scores. Highly synergistic and selective drug target pairs are colored red, while the remaining ones are colored blue. (B) A network of pathway associations 
of the predicted synergistic and highly selective drug target pairs (Supplementary Table S4). Node size represents the number of genes that participate in a pathway. 
Edge width represents the number of synergistic and selective pairs involving genes in the two connected pathways (its nodes). 



potentially damage various liver functions (Supplementary 
Table S3} . Hence, as noted above for single targets, the current 
lack of genome-scale metabolic models for various tissues 
prevents the prediction of other potential sources of cellular 
damage that may be caused by these knockdowns— and it 
would be reasonable to assume that many of these pairs 
would be actually labeled as toxic as new tissue models are 
developed. 



The targeting of a gene whose synthetic lethal 
partner is inactivated in specific cancer types 
leads to selective treatments 

The specific targeting of a gene participating in a synergistic 
pair is especially appealing in tumors in which its interacting 
gene is specifically inactivated. The targeting of such a gene 
solely is likely to selectively damage the tumor, without 
affecting the function of healthy tissues in which the 
interacting gene is not inactivated and hence can still back 
up for the targeted gene. Notably, this strategy can achieve 
therapeutic selectivity without requiring the ability to compu- 
tationally predict the effect of double knockdowns on healthy 
tissues. The inactivation of one of the interacting genes 
specifically in cancer cells may result from cancer type-specific 
regulation of metabohc activity, or from the dysfunction of 
metabohc enzymes due to cancer somatic mutations. Next, we 
utilized genomic and transcriptomic data to infer gene 
inactivation across an array of cancers, which leads to the 
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identification of cancer type-specific targets based on the 
predicted synergistic gene pairs. 

Germline mutations in two metabohc enzymes, succinate 
dehydrogenase (SDH) and fumarate hydratase (FH) have been 
previously found to be causally imphcated in oncogenesis 
(Futreal et aZ, 2004). In both cases, the malignancy was 
characterized by somatic loss of the corresponding wild-type 
allele, resulting in complete loss of enzymatic activity. Our 
model predicts SDH to be synthetic lethal with PC, showing 
that either one of these enzymes is required to satisfy TCA 
anaplerotic demands essential for the synthesis of aspartate 
(that is used to further synthesize arginine and pyrimidine 
nucleotides) from oxaloacetate (OAA) (Figure 4A}. OAAcanbe 
produced either by PC from pyruvate or through glutamino- 
lysis that produces a-ketogluterate (aKG), which is converted 
to OAA through several TCA cycle reactions. SDH dysfunction 
blocks this anaplerotic flow from glutaminolysis toward OAA. 
Hence, inhibition of PC is predicted to specifically inhibit 
proHferation of SDH-deficient cancer cells, without damaging 
normal cells. Several experimental drugs that target PC are 
currently available and may potentially be used for treating 
SDH-deficient cancers such as paragangliomas and pheo- 
chromocytomas (Frezza et aZ, 2011). A similar argument 
regarding synthetic lethality between SDH and PC also holds 
for FH and PC, as FH deficiency is also expected to block 
metabohc flux toward OAA. Our model predicts that FH forms 
synthetic lethality with nine additional genes (Supplementary 
Table S3}, four of which already have drugs that target them. 
These pairs can hence be potentially used to target specific 
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Figure 4 Predicted synthetic lethal gene pairs in TCA cycle (A), phosphatidylserine biosynthesis (B) and pentose-phosphate pathway (C). Red arrows represent 
predicted synthetic lethality. In (A) SDH undergoes loss of heterozygosity resulting in a complete loss-of -function in paragangliomas and pheochromocytomas; 
in (B) PTDSS2, CEPT1 and PCYT2 undergo chromosomal deletions in testicular cancer (PTDSS2), pheochromocytoma (CEPT1) and renal cell carcinoma (PCYT2); 
in (C) RPE undergo chromosomal deletions in several cancer types, including adrenocortical carcinoma and head and neck squamous cell carcinoma. In all cases, 
the targeting of the corresponding synthetic lethal partner is expected to selectively affect proliferation of the cancer cells. 



cancers associated with FH mutations such as leiomyoma, 
leiomyosarcoma or renal cell carcinoma (Futreal et aZ, 2004} 
(though notably, the tumor suppressor function of FH was 
recently shown to go beyond its TCA cycle activity, contri- 
buting to DNA damage response; Yogev et aU 2010). An 
experimental validation of a predicted synergy between 
FH and an enzyme in heme metabohsm has already been 
performed in a follow-up study described in the Discussion. 
Overall, in addition to FH and SDH, somatic mutations were 
previously found in 13 other metabohc enzymes (Forbes et aZ, 
2009). These 13 genes are predicted to participate in 44 
synergistic gene pairs (Supplementary Figure S3} . 

Chromosomal deletions of metabohc genes in specific 
cancer types is expected to lead to either a complete loss of 
function or reduced activity due to dosage effects. Analyzing 
data from the Atlas of Genetics and Cytogenetics in Oncology 
and Haematology (Huret et aU 2004} revealed that chromo- 
somal deletions in all genes that participate in the predicted 
synthetic pairs were previously identified in at least one out of 
103 cancer types. Hence, the targeting of the corresponding 
synthetic lethal genes in these cases is expected to selectively 
damage cancer cells. For example, our model predicts several 
synthetic lethal interactions between enzymes in two alter- 
native pathways for the biosynthesis of phosphatidylserine 
(accounting for 5 to 10% of cell membrane phospholipids}, 
utilizing either ethanolamine or choline as precursors 



(Figure 4B}. These include, synthetic lethality between phos- 
phatidylserine synthase 1 (PTDSSl} (in the choline utilizing 
pathway}, and phosphatidylserine synthase 2 (PTDSS2}, 
ethanolamine phosphotransferase (CEPTl} and ethanol- 
amine-phosphate cytidylyltransferase (PCYT2} (all three in 
the ethanolamine pathway} . The predicted synergy between 
these two pathways was previously demonstrated experimen- 
tally (Arikketh et al, 2008; Vance and Vance, 2009}. We find 
that each of the three enzymes in the ethanolamine pathway 
are frequently deleted in various cancer types (Supplementary 
Table S13}, including testicular cancer (PTDSS2}, pheo- 
chromocytoma (CEPTl} and renal cell carcinoma (PCYT2}, 
suggesting that the targeting of their synthetic lethal partner 
PTDSSl would selectively inhibit prohferation of these specific 
cancers. Another example is that of ribulose 5-phosphate 
3-epimerase (RPE}, participating in the non-oxidative branch 
of the PPP. This enzyme is predicted to be synthetic lethal with 
all three enzymes in the oxidative branch of the pathway, 
glucose 6-phosphate dehydrogenase (G6PD}, 6-phosphoglu- 
conolactonase (PGLS} and phosphogluconate dehydrogenase 
(PGD} (Figure 4C}. The latter is due to the requirement of 
either the oxidative or non-oxidative branches of the PPP 
for synthesizing a-D-ribose 5-phosphate, which is an essential 
precursor for nucleotides. We find that RPE undergoes frequent 
deletions in several cancer types, including adrenocortical 
carcinoma and head and neck squamous cell carcinoma. 
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2 Oral sec 

g Hepatocellular carcinoma 
Colorectal carcinoma 
Spermatocytic seminoma 
Thyroid carcinoma, follicular 
Thymic carcinoma 
Prostate cancer 
Fallopian tube carcinoma 
Esophageal carcinoma 
Gastric cancer 
Testis 

Head and neck SCC 
Uterine/cervix carcinoma 
Osteosarcoma 
Cervical carcinoma 
Adrenocortical carcinoma 
Pheochromocytoma 



Deleted genes 




Drugs 



Figure 5 Drugs predicted to target specific cancer types based on chromosomal loss of synthetic lethal participant genes. Specific drugs may be effective in treating 
specific types of cancer, based on our predictions of selective synthetic lethal gene pairs. Cancer types that show a high frequency (in yellow and white) of chromosomal 
deletions of specific genes are susceptible to drugs inhibiting the genes' synthetic lethal complements. Experimental drugs are followed by an asterisk. In cases where 
multiple drugs target the same gene, only a single representative drug is shown here with the remaining drugs specified in Supplementary Table S9. 



suggesting that the targeting either of its corresponding 
synthetic lethal partners should be useful for selectively 
inhibiting prohferation of these specific cancers. Notably, the 
experimental drug hydroxyacetic acid was previously shown 
to inhibit G6PD activity, suggesting its potential usage for 
treating the above-mentioned cancer types. Figure 5 shows the 
association of approved and experimental drugs with cancer 
types, in accordance with cancer-type chromosomal deletions 
of synthetic lethal genes (Supplementary Table S9} . 

Analyzing gene expression data of healthy human tissues 
taken from Su et al (2004), we have found that out of the 342 
predicted synthetic lethal targets, 72 pairs have at least one 
gene with a significantly low expression level in a target cancer 
type and not in healthy tissues, spanning leukemia, melanoma, 
lung, colon, CNS, ovarian, renal, prostate and breast cancer 
(see Materials and methods). Taken together, these results 
testify to the promising utility of this approach for identifying 
new potential treatment targets for specific cancers. 

Discussion 

As a further experimental support for the plausibility of the 
predicted targets, a follow-up work by Frezza (under review) 
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studied the predicted synergy (described above; Figure 2) 
between the TCA cycle enzyme FH and the heme metabolism 
pathway. This synergy is of a particular interest both from a 
practical perspective, as FH loss-of-function is associated with 
Hereditary Leiomyomatosis and Renal-Cell Cancer (HLRCC) 
(Tomlinson et aZ, 2002) and from a basic science perspective 
as, currently, no mechanism that explains the capability of 
these cells to survive without FH (and hence without a func- 
tional TCA cycle) has been provided. Notably, the synergy with 
heme metabohsm predicted by the cancer model is not evident 
based on the gene expression data measured under FH- 
deficient cells solely, and cannot be gleaned from analyzing 
the human metabohc model of Duarte et al (2007) . In the study 
by Frezza (under review), newly characterized genetically 
modified kidney mouse cells carrying a conditionally targeted 
FH allele were used to confirm that targeting heme metabohsm 
would render FH-deficient cells non-viable while not dama- 
ging the control wild-type FH-containing cells. This was found 
to result from a linear pathway involving the biosynthesis and 
degradation of heme, which enables FH-deficient cells to 
utilize the accumulated TCA cycle metabolites and permits 
partial mitochondrial NADH production. These results confirm 
the suggested targeting of heme metabolism to treat HLRCC 
patients, providing an interesting specific demonstration of 
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the utility of using our generic cancer model, and the usage of 
synthetic lethality as a mean for specifically targeting cancer 
cells. Future validation of predicted targets could involve 
specific shRNA silencing of predicted single and synergistic 
targets, across different cancer cell lines. Additionally, 
predicted repurposing and combinations of approved 
and experimental drugs can be easily tested across cell lines, 
measuring the reduction in prohferation rate following their 
treatment. Finally, of particular interest in this respect are 
'hubs' participating in many predicted SLs, whose knockdown 
could have large functional effects. 

Beyond the prediction of new potential drug treatments, the 
modeling approach presented here is expected to open up 
many additional exciting possibilities in the near future. First, 
data on gene regulation may be integrated with the stoichio- 
metric model using existing approaches (Covert et aU 2004; 
Shlomi et aU 2007} to further enhance the model's accuracy. 
Second, the forthcoming introduction of a new version of the 
human metabohc model may further increase the accuracy 
of the resulting cancer models. Third, the reconstruction of 
human tissue-specific models, such as the recently developed 
models for liver (Jerby et aU 2010} and kidney (Chang et aU 
2010}, can lead to better computational predictors of drug 
selectivity. A major current difficulty in the reconstruction of 
such models involves the definition of a metabohc 'objective 
function' for normal tissues. Fourth, future studies could probe 
the functional effects of the gradual knockdown of reactions, 
aiming to capture a richer repertoire of enzymatic inhibition 
and gene knockdown. Finally, individual cancer cell line 
models may lead to the prediction of specific cancer 
biomarkers, building on existing CBM methods (Shlomi et aU 
2008}. In summary, the model presented here lays down a 
fundamental computational counterpart for interpreting the 
rapidly accumulating proteomics (Bichsel et aU 2001} and 
metabolomics (Fan et aU 2009} data characterizing cancer 
metabohc alterations, and paves the way both for obtaining 
a systems level understanding of cancer metabohsm and for 
designing new therapeutic means that selectively target them. 

Materials and methods 

Reconstructing a human cancer metabolic model 

To reconstruct a cancer metabolic model in humans, we collected 
gene expression data from all cancer cell lines in the NCI-60 collection 
(Grever etal, 1992}, a collection of microarray experiments in 59 cancer 
cell lines of various types performed in a controlled and consistent 
manner with a single platform and within a single laboratory, and 
assembled a core set of 197 cancer genes that are highly expressed (with 
intensity above 200} in >90% of cell line measurements (while the 
results of our analysis remain highly robust to different definitions of 
core gene sets; Supplementary Information}. Assuming that reactions 
associated with these highly expressed cancer genes are likely to be 
activated in cancer cells, we applied the following approach: our goal is 
to extract from the human metabolic model of Duarte et al (2007}, the 
most parsimonious, consistent model, that includes this core reactions 
set and enables cellular proliferation. A metabolic network model is 
considered consistent if it enables the activation all of its reactions — 
that is for each reaction there exist a feasible flux distribution in 
the model (satisfying stoichiometric mass balance and directionality 
constraints} in which this reaction is carrying a non-zero flux. Bearing 
this definition, we employ a greedy search heuristic approach that is 
based on iteratively pruning reactions from the human metabolic 
model in a random order, while maintaining the consistency of the 
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pruned model with regard to its core reactions set (starting from a 
generic version of the human metabolic model that is consistent} . In 
each pruning step, a reaction is removed from the cancer model being 
built only if its removal does not prevent the activation of the reactions 
in the core reactions set. Since the reactions' scanning order may affect 
the resulting model, the algorithm is executed with different, random 
pruning orders (1000 in our implementation} to construct multiple 
candidate models. The fraction of models that contain a certain reaction 
reflects the latter's confidence score, which indicates the extent to 
which it should be included in the final cancer model. Hence, to 
construct the final cancer model, the candidate models are aggregated, 
starting from core reactions set and iteratively adding reactions ordered 
by their confidence scores, until a final, minimal but consistent with the 
core reactions, model is obtained. A similar approach has been 
concurrently used to develop a model of healthy liver metabolism 
(Jerby et al, 2010} (Supplementary Information}. Notably, blocked 
reactions from the model of Duarte et al were removed before the 
application of the MBA algorithm. 

To predict gene contribution to biomass production, we added a 
growth reaction to the resulting model, representing the steady-state 
consumption of biomass compounds required for cellular proliferation. 
The stoichiometric coefficients of the growth reaction represent the 
relative molecular concentrations of 42 essential metabolites, including 
nucleotides, deoxynucleotides, amino acids, lipids, etc. in human 
tissues. These relative concentrations are calculated based on data 
regarding mass composition of liver and muscle tissues (Supplemen- 
tary Table S6}. A sensitivity analysis shows that the prediction 
performance of the cancer model is highly insensitive to the specific 
definition of biomass composition (Supplementary Information}. In all 
simulations, we assume a standard RPMI-1640 medium, in accordance 
with the medium used in our reference shRNA experimental data set 
(Luo et al, 2008; Supplementary Table S7} . The biomass production rate 
predicted by the cancer model is 40 % lower than that predicted by the 
human network model, reflecting that the two models are indeed 
functionally different. Notably, the generic human model does not 
represent a concrete cell-type (but rather a collection of reactions that 
take place within different cell types} , and hence its predicted optimal 
biomass production rate does not accurately represent a rate that is 
achievable by a specific cell type. By definition, the maximal biomass 
production rate in the cancer model cannot be higher than that 
achievable in the generic human model as the cancer model consists of 
a subset of the reactions of generic model. 

FBA was used to simulate the effects of gene knockdowns on the 
biomass production rates of the proliferating cells (Forster et al, 2003}. 
Genes whose knockdown reduces the maximal biomass production 
rate in > 1 % were considered to be growth supporting. Notably, the 
results presented in the paper are insensitive to the specific choice 
of threshold for defining growth reduction, as the majority of the 
knockouts (99.13%} result in either no reduction of growth or its 
complete elimination (Supplementary Information} . 



Validating predicted growth-supporting genes 
based on shRNA data 

Luo et al (2008} provide a ranking of gene probes that are targeted by 
shRNA, based on the effect that their knockdown exerts on cancer 
cellular proliferation. To check whether a set of predicted growth- 
supporting genes are ranked as highly essential for cancer proli- 
feration, we computed a one-sided KS statistic, comparing the ranking 
distribution of the predicted gene set probes with the ranking 
distribution of the remaining probes. Significance is assessed by 
comparing the derived KS statistic with similar statistics computed for 
10 000 randomly selected gene sets of the same size (in accordance 
with the statistical method underlying Gene Set Enrichment Analysis 
statistic; Subramanian et al, 2005}. 

Computing cytostatic scores for single and double 
drug target predictions 

A cytostatic score of a gene represents the extent to which its 
knockdown obliterates cancer growth compared with its effect on 
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ATP production in the normal, human metabolic model. The gene 
knockdown effect on growth rate is computed by applying FBA on the 
cancer model, denoting by WTgrowth and KOgrowth the growth rate 
before and after the knockdown, respectively. The gene knockdown 
effect on ATP production is computed by applying FBA on the human 
metabolic model, measuring the reduction in maximal achievable 
ATP production rate, denoting by WTatp and KOatp the maximal 
ATP production rate before and after the knockdown, respectively. The 
cytostatic score is defined as (KOatp/WTatp) x (l-KOgj-owth/WTgj-owth)- 
A cytostatic score of 1 denotes a highly selective knockdown that 
completely eliminates cancer growth without affecting ATP production 
of healthy cells, and a score of 0 denotes a highly non-selective 
knockdown. The cytostatic score for synergistic gene pairs is computed 
in an analogous manner, by simulating the effect of double knock- 
downs on growth and ATP production rates. 



Predicting and validating synergistic drug targets 

We define a synergy score for a gene pair as the additional drop in 
growth rate below the minimal rate achieved by the single knock- 
downs. Specifically, denoting by KOa, KOb and KOab, the growth rates 
following the knockout of gene A, gene B and the joint knockout 
of genes A and B, respectively, the synergy score is defined as: KOab/ 
min(KOA, KOb) (in accordance with the commonly used 'highest single 
agent' synergy model; Berenbaum, 1989). The validation of the 
predicted synergistic gene targets via yeast orthologs was based on an 
orthology mapping by Berglund et al (2007). Among the yeast genes 
that participate in the high-throughput genetic interaction screen 
(Costanzo et al, 2010), we identified 75 genes that are orthologous to 
human genes that participate in predicted synergistic targets. 

The correlation between drug efficacy and gene expression was 
computed as following: for each gene pair in which one gene is 
targeted by drug D and the other one is denoted G, we computed the 
Spearman correlation between the effective growth inhibition con- 
centrations (GI-50) of D with the expression levels of G across the 60 
cell lines. The P-value was computed by a Wilcoxon test, comparing 
the distribution of Spearman correlations for all predicted synergetic 
drug-gene pairs, with the distribution of Spearman correlations for a 
large random sample of drug-gene pairs. 

To identify synergistic gene pairs that consist of genes with signi- 
ficantly low expression levels specifically in cancer cells, we obtained 
gene expression data from 75 healthy human tissues (Su etal, 2004) and 
compared it with the NCI-60 cancer cell line expression data (Grever 
et al, 1992). First, we normalized the tissue expression data set to the 
same mean and s.d. as the NCl-60 cancer cell line expression data. Then, 
we applied Wilcoxon test to identify genes with significantly low 
expression level in at least one out of eight cancer types represented in 
NCI-60, followed by false discovery rate correction for multiple testing. 



Supplementary information 

Supplementary information is available at the Molecular Systems 
Biology website (www.nature.com/msb). 
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